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The flow over afterbody geometries has been investigated' using the reduced 
Navier-Stokes (RNS) approximation. Both pressure velocity flux-split and composite 
velocity primitive variable formulations have been considered. Pressure or pseudo- 
potential relaxation procedures are combined with sparse matrix or coupled strongly 
implicit algorithms to form a 3-D solver for general non- orthogonal coordinates. 
Three-Dimensional subsonic and transonic viscous/inviscid interacting flows have 
been evaluated. Solutions with and without regions of recirculation have been 
obtained. 
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INTRODUCTION 


For a significant class of problems, exhibiting strong viscous - inviscid 
interaction, the reduced form of the Navier-Stokes (RNS) equations can provide an 
effective alternative to the full Navier-Stokes (NS) system for large Reynolds (Re) 
number numerical calculations, in both, two- and, three- dimensions. In a series of 
papers [1-4] , Rubin and Khosla have developed and applied composite velocity 
primitive variable and pressure velocity formulations for a variety of sub-, trans- 
and supersonic flow problems. The RNS formulation, which represents a composite of 
Euler, higher-order boundary- layer triple deck equations, is consistent with all 
asymptotic large Re strong interaction theories and results in a uniformly valid 
single set of equations. The RNS system includes all normal and secondary flow 
diffusion terms in the surface momentum equations, but does not include normal or 
axial flow diffusion effects in the surface normal momentum equation. This 
approximation is self consistent for the large Re limit. Crossflow or secondary flow 
diffusion is retained in all equations as these effects are required to accurately 
model many three dimensional external and internal flows; see e.g. references [5- 
11], where it has been shown that explicitly added artificial viscosity can 
significantly destort secondary flow, heat transfer and surface stress predictions. 

In the present paper, both the pressure and composite velocity primitive 
variable RNS formulations are considered for the computation of three-dimensional 

afterbody flows. In the pressure velocity formulation, for the velocity (q) and the 
pressure (p) , the dominant elliptic behavior and associated upstream influence, as 
manifested through an eigenvalue analysis [12] , is associated with acoustic 
propagation and therefore the discretization requires an appropriate form of 
pressure velocity flux-split differencing. In the composite velocity formulation, 
the upstream influence is modelled by considering the velocity as a multiplicative 
composite of viscous or rotational velocities (U,W) , and velocities derived from a 
pseudo potential ($) . This representation of velocities is defined in the spirit of 
matched asymptotic expansions and is such that the resulting system reduces to the 
potential equation for inviscid irrotational regions. For rotational inviscid flows, 
a composite form of the Euler equations is recovered and a vorticity (U,W) term 
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appears on the right hand side of the pseudo-potential continuity equation. The 
composite formulation compliments the two-dimensional vorticity stream function 
formulation, but is more suitable for three -dimens ion problems, as the total number 
of unknowns is does not increase, as would be the case for the three-dimensional 
vorticity- stream function procedure. 

Flux vector splitting [12] is the basis of the discretization for the (q,p) 
system and flux biasing with appropriate upwinding [1] is specified for the (U,W,$) 
system. For supersonic flows the both RNS systems reduce to initial value PNS 

formulations. The (q,p) system can be solved with a spatial "boundary- layer" type 
marching relaxation procedure. Since flux biasing is employed for the composite 
(U,W,$) system, supersonic flows are computed by a relaxation process with the 
Enquist-Osher compressibilty correction applied for the pseudo-potential $. It 
should be emphasized that both formulations are valid throughout the Mach number 
Reynolds number range, including M £o « 1 and M » 1, and do not require the addition 

of any explicitly added artificial viscosity. The inherent numerical viscosity 
associated with the discretization is sufficient to capture strong shocks over 3 to 
4 mesh points. This numerical error is minimized with fine grids and a multi-grid 
strategy [1]. In earlier investigations of Rubin and Khosla, both pressure and 
composite velocity RNS codes have been applied for the computation of 
incompressible, transonic and supersonic flows with strong viscous- inviscid 
interaction, flow reversal and shock capturing e.g. see references [1-11]. 
Furthermore, the omitted diffusion terms have been incorporated directly, or via a 
deferred corrector approach, to obtain solutions of the complete NS equations. It 
has been shown that these effects are minimal for the problems considered. These 
procedures have been tested for the steady flow over an axisymmetric boattail and 
for the unsteady flow over a Joukowski airfoil. In the steady RNS formulations, 
outflow boundary conditions are required only for the pressure or pseudo potential 
and, as such, these are the only variables for which global storage of velocities is 
required outside of reverse flow regions. This is a significant simplification for 
steady three dimensional flow computations, and this can also be useful in the 
computation of complex three-dimensional flows on very fine meshes. 

In both RNS formulations, the axial convective terms are first or second- order 
upwind differenced. The solution procedure takes advantage of the flux split/upwind 
differencing and results in a boundary -layer type streamwise marching method that is 
imbedded in a global relaxation process. In certain strong interaction cases a 
sparse matrix direct solver is applied for the pressure variable cross plan 
solution; in the composite velocity velocity the coupled strongly implicit ILU 
inversion is applied. 

Subsonic flow past afterbody configurations with elliptic and hyperelliptic 
cross sections are discussed in this paper. Comparison are given between the 
pressure and composite velocity solutions. Grid resolutions studies are used to 
assess the accuracy of the two solution procedures. 


GOVERNING RNS SYSTEM AND SOLUTION PROCEDURE 

As described in earlier investigations [1-11], the RNS approximation leads to a 
single composite system that includes the Euler, second-order boundary layer and 
triple deck equations. The governing RNS system is considered for three-dimensional 
generalized coordinates and is designed to allow for shock capturing and flow 
reversal. In nondimensionalized form, the RNS system for low- speed flow is given as 
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Continuity Equation 


~ (pjg u) + “ (pj g v) + ~ (pjg w) = 0 (1-a) 

£ - Momentum Equation 

(pjg u ) + (pjg VU) + J^~(pjg wu) + curvature terms = 

- g 11 p^ - g 12 p^ - g 13 p^. + viscous terms (1-b) 

f - Momentum Equation 

d 3 ^ 2 

— (pjg wu) + — (pTi vw) + J^~(pJ g w ) + curvature terms = 

' g 3 1 " g 32 • g 33 P^. + viscous terms (1-c) 


tj - Momentum Equation 

3 3 3 

— (pjg uv) + (pjg V 2 ) + (pjg vw) + curvature terms = 
- g 21 P^ * g 22 P^ - g 23 P f 
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where g 1 ^ are the contravariant form of the metrics, Jg is the Jacobian , 
p,u,v,w,p and T are the density, velocities, pressure and temperature, repectively 
and 7 is the ratio of specific heats. In these equations all variables are non- 
dimensionalized with respect to their free stream conditions except for the pressure 
which is normalized with the free stream dynamic pressure. 


a) Pressure/Velocity Formulation (P.V.) 
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The pressure flux vector splitting, as described in reference [12], leads to the 
following discrete representation of the axial pressure gradient: 

V w i-l/2 (p i ' P i-1 } 7 + (lw i+l/2 ) (p i+l ' P i )/Z ^ (2) 

where 

w < = min [7 M /(l+( 7 -l)M|) , 1 ] for constant stagnation enthalpy 

and 


w < = min[M 2 , 1] for the full flux split energy equation. 


This pressure gradient splitting and the associated convective upwinding satisfies 
the major eigenvalue continuity constraints on the fluxes and flux derivatives. For 
w = one eigenvalue is alwaya zero so that sharp shocks are obtained. For regions 

of reverse, flow the condition w=0 is required. This ensures that the fluxes, flux 
derivatives and eigenvalues remain continuous throughout the flow. This form of flux 
splitting is designed to maintain a bias in the direction of convective fluxes and 
leads to a relaxation procedure that is solely acoustically driven in the subsonic 
regions . 

As described in earlier investigations, the continuity equation is discretized 
at (i,j-l/2,k). The streamwise (£- ) and crossflow (f- ) momentum equations are 

discretized at (i,j,k) and the normal (ij- ) momentum equation is discretized at 
(i,j+l/2,k). This discretization is consistent with the flux eigenvalues and with 
the appropriate boundary conditions required for u and p . For outflow boundaries 
without recirculation, only a pressure or pressure gradient condition is required. 
Far field boundary conditions, with positive outflow, are only required for u, p, 
and T. Inflow boundary conditions are specified for all flow variables or their 
gradients. Zero injection is assumed for all solid surface. The resulting algebraic 

system for the delta form of pressure Ap=p n+ ^- p n is solved in each cross-plane using 
a coupled version of a sparse matrix direct solver. Additional details on the solver 
can be found in references [11,14] 

b) Composite Velocity (U,W,<!>) Formulation (C.V.) 

In the spirit of matched asymptotic expansions, the contravariant velocities are 
rewritten as : 

u = (U+l) ( g + g ^ + g = (U+l)u e 

21 22 23 

V = (g + g ^ + g (3) 

13 32 33 

w - (U+l)(g + g fc + g $ ) + W 

The composite representations for u and w, the axial and cross -flow velocity 
components, contain two types of terms, e.g. a rotational "pseudo" potential 
function <E> and viscous velocities U, W. The viscous no-slip boundary conditions are 
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introduced through the velocities U, W. The kinematical boundary conditions are 
satisfied through the pseudo potential. This is consistent with the potential and 
boundary- layer approximations. Substitution of equation (3) into equation (1) leads 
to the RNS composite system. Additional details and the resulting equations are 
given in references [1-4] . All derivatives are approximated using three point 
central differences except for the convective U^ , terms: these are second- order 

upwinded. Upstream influence in attached flow regions arises solely through $ and 

thus an outflow boundary condition on <3?, in lieu of p for the pressure formulation, 
is required. This leads to boundary- layer like streamwise marching for U, W with 
global relaxation for $ . In reversed flow regions , w^ is upwind differenced and 

therefore the relaxation procedure includes both 3> and w, as is the case for p and 
(u,w) in the pressure split velocity formulation. The resulting algebraic system in 


c n+l 




each crossplane, for the delta form of $, V<$=$ - $ is solved using a consistent 
version of the coupled strongly implicit algorithm. This algorithm has previously 
been described in reference 10. 


BOUNDARY CONDITIONS 

The boundary conditions used for this investigation are; (i) uniform flow at the 
inflow; (ii) weak viscous/inviscid interaction at the outflow, thus at £ = 

0 or the negative eigenvalue fluxes are set to zero for the pressure velocity 
solver; (iii) the no slip and zero injection conditions are specified at the body 
surface; i.e. u = U+l =0, v = v g = 0 and w = W = 0; (iv) far from the body the flow 

is assumed to be undisturbed; i.e. u = 1, w=0, $>= , T = 1 and the entropy is 

assumed to be constant. 


RESULTS 

Afterbody configurations of elliptical (Figure la) and hyper-elliptical (Figure 
lb) cross-sections have been examined. In view of the geometric symmetry, the 
calculations have been carried out in one quarter of the flowfield. A body fitted 
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Figure 1-a: Surface grid distribution 
- ellipse 



Figure 1-b: Surface grid distribution 
- hyperellipse 
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computational grid e.g., see Figure 2 for an elliptic cross-section, was generated 


Figure 2 : 

Ellipse - grid in the cross -plane 71x17 

An = 0.00175 
mm 

variable stretching 


T 1 

5 6 



by using a shearing transformation. A stretched grid is employed in the normal 
direction while a uniform spacing is prescribed in the azimuthal direction. The 
axial boattail is generated by fitting a cubic polynomial between two appropriate 
cross-sections . 

Both laminar and turbulent flows have been investigated using the composite 
velocity (C.V.) formulation. Only laminar flow computations have been carried out 
using the pressure velocity (P.V.) formulation. Figure 3a and 3b depict the pressure 




Figure 3 -a: Pressure Coefficient 


Figure 3-b: Skin-Friction Parameter 


- comparison P.V. vs C.V. 


-comparison P.V. vs C.V. 


coefficient and skin friction parameter for a laminar non- separated boattail of 
elliptical cross-section using both (P.V.) and (C.V.) formulations. The resonable 
agreement between the two solutions provides a self consistent evaluation of the two 
formulations. The metrics are calculated analytically for the P.V. solutions and 
numerically for the C.V. solutions. As the grid is refined the two solutions merge. 
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The pressure coefficient and skin friction parameter for turbulent flow 
computations for an afterbody of elliptical cross-section at M =0 . 9 using the (C.V.) 

formulation are shown in Figures 4-a and 4-b. As seen in the figure, a sharp weak 
transonic shock has been captured quite well. 

The crossplane grid for an afterbody of hyper-ellipse cross-section is depicted 

4 

in Figure 5. A comparison of the laminar P.V. flow solution for Re=10 , on two grids 



x 

Figure 4-a: Pressure Coefficient 
- ellipse 



- ellipse 



Figure 5: Hyperellipse 
grid in the cross-plane 79x17 


A ' , uiin“° ' 0005 


Stretching factor a=1.09 


61x79x17 and 121x79x17, is shown in Figures 6a and 6b. The pressure coefficient and 
skin friction parameter along a symmetry plane and the location of maximum azimuthal 
curvature are used for this comparison. A region of large recirculation, along the 
boattail, has been computed by this technique. The sensitivity to grid refinement is 
quite evident. Further grid refinement in both the axial and azimuthal directions 
are most likely required to obtain an acceptable level of accuracy. The axial 
velocity vectors in the recirculation region are depicted in Figures 7a and 7b. The 
effect of different subsonic Mach numbers is shown in Figures 8a and 8b. The 
recirculation is slightly increased. This is consistent with previous subsonic 
calculation on axisymmetric boattails [2] . 
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Pressure Coefficient 












From the above disscusion, it is clear that the present RNS techniques have the 
ability to compute fairly complex flow fields. The computational cost for the C.V. 
solution is modest, e.g., with 81923 grid points and 5 coupled unknowns this is of 
the order of 30 to 60 minutes on the Cray Y-MP. The computational cost of the P.V. 
solutions, on the other hand, is much greater. This is due of the fact that that 
sparse matrix direct solver, which requires 8 to 9 times the CPU of the CSIP, is 
applied for much of the cross plane inversion. Strategies to reduce the overall 
computational effort by minimizing the need for multiple LU (direct solver) 
inversions have been investigated. The computational times can be reduced 
considerably with bach subsitution. Recent studies have provided factors of less 
than 2 over C.V. computer time . A stabilized version of the consitent coupled 
strongly implicit algorithm for the P.V. formulation is also currently under 
investigation. This will result in further reductions in the P.V. CPU time. 


CONCLUSION 

Two formulations of the three-dimensional RNS equations have been investigated 
for the computation of laminar/turbulent subsonic , transonic flows with large 
recirculation regions. The technques are quite efficient in terms of storage and 
computational times. Future applicaltions to complex internal/external geometries 
and supersonic free streams are under investigation. 


ACKNOWLEDGEMENT 

The work has been supported by the AFOSR (L. Sakell, Technical Monitor) under 
Grant No. AFOSR 90-0096. 

The computations have been performed on the Cray Y-MP through a project 
supported by the Ohio Supercomputer Center. 

REFERENCES 

1. Gordnier, R.E. and Rubin, S.G., "3-D Composite Velocity Solutions for 
Subsonic/Transonic Flows over Afterbodies", AIAA paper no. 89-1837, Presented at the 
20th Fluid Dynamics, Plasma Dynamics and Lasers Conference, Buffalo, New York, June 
12-14, 1989. 

2. Swanson, R.C., Rubin, S.G. and Khosla, P.K. , "Calculation of Afterbody Flows 
with a Composite Velocity Formulation", AIAA paper no. 83-1736, Presented at the 
AIAA 16th Fluid and Plasma Dynamics Conference, Danvers, Massachusetts, July 12- 
14,1983. 

3. Khosla, P.K. and Rubin, S.G., "A Composite Velocity Procedure for the 
Compressible Navier-Stokes Equations", AIAA Journal, Vol. 21, No. 11, November 1983, 
pp. 1546-1551. 

4. Rubin, S.G., Celestina, M. and Khosla, P.K. , "Second Order Composite Velocity 
Solution for Large Reynolds Number Flows", AIAA paper no. 84-0172, Presented at the 
AIAA 22nd Aerospace Sciences Meeting, Reno, Nevada, January 9-12, 1984. 

5. Reddy, D.R. and Rubin, S.G., "Consistent Boundary Conditions for Reduced 
Navier-Stokes (RNS) Scheme Applied to Three-Dimensional Viscous Flows", Journal of 
Fluids Engineering, Vol. 110, September 1988, pp. 306-314.. 


591 



6. Rubin, S.G. and Reddy, D.R., "Analysis of Global Pressure Relaxation for 
Flows with Strong Interaction and Separation", Computers and Fluids, Vol. 11, 1983, 
pp . 281-306. 

7. Khosla, P.K. and Lai, H.T. , "Global PNS Solutions for Transonic Strong 
Interaction Flows", AIAA 22nd Aerospace Science Meeting, AIAA paper no. 84-0458, 
January, 1984. 

8. Cohen, R. and Khosla, P.K., "Three-Dimensional Reduced Navier Stokes 
Equations solution for Subsoinc Separated and Non- Separated Flows Using a Global 
Pressure Relaxation Procedure", International Journal for Numerical Methods in 
Fluids, Vol. 9, 1989, pp . 1087-1098. 

9. Rubin, S.G. and Khosla, P.K., "Navier -Stokes Calculations with a Coupled 
Strongly Implicit Procedure", Comput. Fluids 9, 1979, pp. 163-180. 

10. Khosla, P.K. and Rubin, S.G., "Consistent Strongly Implicit Iterative 
Procedures for Two-Dimensional Unsteady and Three-Dimensional Space -Marching Flow 
Calculations", Computers and Fluids, Vol. 15, No. 4, 1987, pp. 361-377. 

11. Bender, E.E. and Khosla, P.K. , "Two-Dimensional Navier-Stokes Equations using 
Sparse Matrix Solvers", AIAA paper no. 87-0603, Presented at the 25th Aerospace 
Sciences Meeting, Reno, Nevada, January 12-15, 1987. 

12. Rubin, S.G., "RNS/Euler Pressure Relaxation and Flux Vector Splitting", 
Technical Note, Computers & Fluids, Vol. 16, No. 4, 1988, pp . 485-490. 

13. Pordal, H.S., Khosla, P.K. and Rubin, S.G., "A Flux-Split Solution Procedure 
for Unsteady Inlet Flows", AIAA paper no. 90-0585, Presented at the 28th Aerospace 
Sciences Meeting, Reno, Nevada, January 8-11. 1990. 

14. Himansu, A., Khosla, P.K. and Rubin, S.G., "Relaxation/Sparse -Matrix Solvers 
for Three-Dimensional Recirculating Flows", AIAA paper no. 89-0552, Presented at the 
27th Aerospace Sciences Meeting, Reno, Nevada, January 9-12, 1989. 


592 



